Para comprender el concepto de correlación, vamos a analizar los datos disponibles sobre la pandemia del Covid-19 a nivel internacional. Para ello, utilizaremos los datos que elabora el proyecto Our World in Data, de la Universidad de Oxford (Reino Unido). Si haces click en el link anterior, podrás entrar al sitio en inglés, pero también descargar los datos y trabajarlos directamente desde el tercer link Download Dataset. Para ahorrar tiempo, puedes ejecutar la siguiente sintaxis y obtener los datos directamente desde su repositorio de GitHub, con datos actualizados diariamente (ATENCIÓN: la actualización implica que los datos pueden presentar variaciones de un momento a otro):
covid <- read.csv(file = ("https://raw.githubusercontent.com/owid/covid-19-data/master/public/data/owid-covid-data.csv"), header=TRUE, sep=",", fileEncoding = "UTF-8")
NOTA: los datos incluídos acá son los reportados hasta el día 22 de marzo a las 10:00 hrs. (Chile), los que son actualizados constantemente
El primer paso para analizar cualquier tipo de datos estadísticos, es utilizar las herramientas de la Estadística Descriptiva, que usted ya conoce. Pero antes, es importante conocer más detalles sobre el objeto de R sobre el cual realizaremos el análisis.
# a) Tipo de objeto (pueden hacerlo para cada uno de los objetos que vayan creando)
class(covid)
## [1] "data.frame"
# b)Variables que contiene (sólo para data.frame)
names(covid)
## [1] "iso_code"
## [2] "continent"
## [3] "location"
## [4] "date"
## [5] "total_cases"
## [6] "new_cases"
## [7] "new_cases_smoothed"
## [8] "total_deaths"
## [9] "new_deaths"
## [10] "new_deaths_smoothed"
## [11] "total_cases_per_million"
## [12] "new_cases_per_million"
## [13] "new_cases_smoothed_per_million"
## [14] "total_deaths_per_million"
## [15] "new_deaths_per_million"
## [16] "new_deaths_smoothed_per_million"
## [17] "reproduction_rate"
## [18] "icu_patients"
## [19] "icu_patients_per_million"
## [20] "hosp_patients"
## [21] "hosp_patients_per_million"
## [22] "weekly_icu_admissions"
## [23] "weekly_icu_admissions_per_million"
## [24] "weekly_hosp_admissions"
## [25] "weekly_hosp_admissions_per_million"
## [26] "new_tests"
## [27] "total_tests"
## [28] "total_tests_per_thousand"
## [29] "new_tests_per_thousand"
## [30] "new_tests_smoothed"
## [31] "new_tests_smoothed_per_thousand"
## [32] "positive_rate"
## [33] "tests_per_case"
## [34] "tests_units"
## [35] "total_vaccinations"
## [36] "people_vaccinated"
## [37] "people_fully_vaccinated"
## [38] "new_vaccinations"
## [39] "new_vaccinations_smoothed"
## [40] "total_vaccinations_per_hundred"
## [41] "people_vaccinated_per_hundred"
## [42] "people_fully_vaccinated_per_hundred"
## [43] "new_vaccinations_smoothed_per_million"
## [44] "stringency_index"
## [45] "population"
## [46] "population_density"
## [47] "median_age"
## [48] "aged_65_older"
## [49] "aged_70_older"
## [50] "gdp_per_capita"
## [51] "extreme_poverty"
## [52] "cardiovasc_death_rate"
## [53] "diabetes_prevalence"
## [54] "female_smokers"
## [55] "male_smokers"
## [56] "handwashing_facilities"
## [57] "hospital_beds_per_thousand"
## [58] "life_expectancy"
## [59] "human_development_index"
# c) Cantidad de columnas (variables) (sólo para data.frame)
ncol(covid)
## [1] 59
# d) Cantidad de filas (casos) (sólo para data.frame)
nrow(covid)
## [1] 77954
# e) Cantidad de filas y columnas (sólo para data.frame)
dim(covid)
## [1] 77954 59
# f) Información sobre las variables (sólo para data.frame)
str(covid)
## 'data.frame': 77954 obs. of 59 variables:
## $ iso_code : chr "AFG" "AFG" "AFG" "AFG" ...
## $ continent : chr "Asia" "Asia" "Asia" "Asia" ...
## $ location : chr "Afghanistan" "Afghanistan" "Afghanistan" "Afghanistan" ...
## $ date : chr "2020-02-24" "2020-02-25" "2020-02-26" "2020-02-27" ...
## $ total_cases : num 1 1 1 1 1 1 1 1 2 4 ...
## $ new_cases : num 1 0 0 0 0 0 0 0 1 2 ...
## $ new_cases_smoothed : num NA NA NA NA NA 0.143 0.143 0 0.143 0.429 ...
## $ total_deaths : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_deaths : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_deaths_smoothed : num NA NA NA NA NA 0 0 0 0 0 ...
## $ total_cases_per_million : num 0.026 0.026 0.026 0.026 0.026 0.026 0.026 0.026 0.051 0.103 ...
## $ new_cases_per_million : num 0.026 0 0 0 0 0 0 0 0.026 0.051 ...
## $ new_cases_smoothed_per_million : num NA NA NA NA NA 0.004 0.004 0 0.004 0.011 ...
## $ total_deaths_per_million : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_deaths_per_million : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_deaths_smoothed_per_million : num NA NA NA NA NA 0 0 0 0 0 ...
## $ reproduction_rate : num NA NA NA NA NA NA NA NA NA NA ...
## $ icu_patients : num NA NA NA NA NA NA NA NA NA NA ...
## $ icu_patients_per_million : num NA NA NA NA NA NA NA NA NA NA ...
## $ hosp_patients : num NA NA NA NA NA NA NA NA NA NA ...
## $ hosp_patients_per_million : num NA NA NA NA NA NA NA NA NA NA ...
## $ weekly_icu_admissions : num NA NA NA NA NA NA NA NA NA NA ...
## $ weekly_icu_admissions_per_million : num NA NA NA NA NA NA NA NA NA NA ...
## $ weekly_hosp_admissions : num NA NA NA NA NA NA NA NA NA NA ...
## $ weekly_hosp_admissions_per_million : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_tests : num NA NA NA NA NA NA NA NA NA NA ...
## $ total_tests : num NA NA NA NA NA NA NA NA NA NA ...
## $ total_tests_per_thousand : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_tests_per_thousand : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_tests_smoothed : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_tests_smoothed_per_thousand : num NA NA NA NA NA NA NA NA NA NA ...
## $ positive_rate : num NA NA NA NA NA NA NA NA NA NA ...
## $ tests_per_case : num NA NA NA NA NA NA NA NA NA NA ...
## $ tests_units : chr "" "" "" "" ...
## $ total_vaccinations : num NA NA NA NA NA NA NA NA NA NA ...
## $ people_vaccinated : num NA NA NA NA NA NA NA NA NA NA ...
## $ people_fully_vaccinated : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_vaccinations : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_vaccinations_smoothed : num NA NA NA NA NA NA NA NA NA NA ...
## $ total_vaccinations_per_hundred : num NA NA NA NA NA NA NA NA NA NA ...
## $ people_vaccinated_per_hundred : num NA NA NA NA NA NA NA NA NA NA ...
## $ people_fully_vaccinated_per_hundred : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_vaccinations_smoothed_per_million: num NA NA NA NA NA NA NA NA NA NA ...
## $ stringency_index : num 8.33 8.33 8.33 8.33 8.33 ...
## $ population : num 38928341 38928341 38928341 38928341 38928341 ...
## $ population_density : num 54.4 54.4 54.4 54.4 54.4 ...
## $ median_age : num 18.6 18.6 18.6 18.6 18.6 18.6 18.6 18.6 18.6 18.6 ...
## $ aged_65_older : num 2.58 2.58 2.58 2.58 2.58 ...
## $ aged_70_older : num 1.34 1.34 1.34 1.34 1.34 ...
## $ gdp_per_capita : num 1804 1804 1804 1804 1804 ...
## $ extreme_poverty : num NA NA NA NA NA NA NA NA NA NA ...
## $ cardiovasc_death_rate : num 597 597 597 597 597 ...
## $ diabetes_prevalence : num 9.59 9.59 9.59 9.59 9.59 9.59 9.59 9.59 9.59 9.59 ...
## $ female_smokers : num NA NA NA NA NA NA NA NA NA NA ...
## $ male_smokers : num NA NA NA NA NA NA NA NA NA NA ...
## $ handwashing_facilities : num 37.7 37.7 37.7 37.7 37.7 ...
## $ hospital_beds_per_thousand : num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...
## $ life_expectancy : num 64.8 64.8 64.8 64.8 64.8 ...
## $ human_development_index : num 0.511 0.511 0.511 0.511 0.511 0.511 0.511 0.511 0.511 0.511 ...
PREGUNTAS a) ¿Cuántas variables contiene la base de datos? b) ¿Cuántos casos o filas contiene la base de datos? c) ¿Cuál es la unidad de análisis?
Considerando que tenemos una base de datos que contiene mucha información, podemos utilizar el paquete dplyr() para filtrar la base de datos sólo a los países de Sudamérica (NOTA: se recomienda fuertemente revisar la nueva base datos con los comandos de la sección anterior, para verificar que se haya filtrado correctamente).
library(dplyr)
sud.covid <- filter(covid, continent == "South America")
sud.covid <- filter(sud.covid, location == "Argentina" | location == "Bolivia" | location == "Brazil" | location == "Chile" | location == "Colombia" | location == "Ecuador" | location == "Paraguay" | location == "Peru" | location == "Uruguay" | location == "Venezuela")
str(sud.covid) #Para conocer la estructura del nuevo conjunto de datos
## 'data.frame': 4003 obs. of 59 variables:
## $ iso_code : chr "ARG" "ARG" "ARG" "ARG" ...
## $ continent : chr "South America" "South America" "South America" "South America" ...
## $ location : chr "Argentina" "Argentina" "Argentina" "Argentina" ...
## $ date : chr "2020-01-01" "2020-01-02" "2020-01-03" "2020-01-04" ...
## $ total_cases : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_cases : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_cases_smoothed : num NA NA NA NA NA NA NA NA NA NA ...
## $ total_deaths : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_deaths : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_deaths_smoothed : num NA NA NA NA NA NA NA NA NA NA ...
## $ total_cases_per_million : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_cases_per_million : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_cases_smoothed_per_million : num NA NA NA NA NA NA NA NA NA NA ...
## $ total_deaths_per_million : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_deaths_per_million : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_deaths_smoothed_per_million : num NA NA NA NA NA NA NA NA NA NA ...
## $ reproduction_rate : num NA NA NA NA NA NA NA NA NA NA ...
## $ icu_patients : num NA NA NA NA NA NA NA NA NA NA ...
## $ icu_patients_per_million : num NA NA NA NA NA NA NA NA NA NA ...
## $ hosp_patients : num NA NA NA NA NA NA NA NA NA NA ...
## $ hosp_patients_per_million : num NA NA NA NA NA NA NA NA NA NA ...
## $ weekly_icu_admissions : num NA NA NA NA NA NA NA NA NA NA ...
## $ weekly_icu_admissions_per_million : num NA NA NA NA NA NA NA NA NA NA ...
## $ weekly_hosp_admissions : num NA NA NA NA NA NA NA NA NA NA ...
## $ weekly_hosp_admissions_per_million : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_tests : num 4 43 4 17 10 53 28 5 1 76 ...
## $ total_tests : num 4 47 51 68 78 131 159 164 165 241 ...
## $ total_tests_per_thousand : num 0 0.001 0.001 0.002 0.002 0.003 0.004 0.004 0.004 0.005 ...
## $ new_tests_per_thousand : num 0 0.001 0 0 0 0.001 0.001 0 0 0.002 ...
## $ new_tests_smoothed : num NA NA NA NA NA NA NA 23 17 27 ...
## $ new_tests_smoothed_per_thousand : num NA NA NA NA NA NA NA 0.001 0 0.001 ...
## $ positive_rate : num NA NA NA NA NA NA NA NA NA NA ...
## $ tests_per_case : num NA NA NA NA NA NA NA NA NA NA ...
## $ tests_units : chr "people tested" "people tested" "people tested" "people tested" ...
## $ total_vaccinations : num NA NA NA NA NA NA NA NA NA NA ...
## $ people_vaccinated : num NA NA NA NA NA NA NA NA NA NA ...
## $ people_fully_vaccinated : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_vaccinations : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_vaccinations_smoothed : num NA NA NA NA NA NA NA NA NA NA ...
## $ total_vaccinations_per_hundred : num NA NA NA NA NA NA NA NA NA NA ...
## $ people_vaccinated_per_hundred : num NA NA NA NA NA NA NA NA NA NA ...
## $ people_fully_vaccinated_per_hundred : num NA NA NA NA NA NA NA NA NA NA ...
## $ new_vaccinations_smoothed_per_million: num NA NA NA NA NA NA NA NA NA NA ...
## $ stringency_index : num 0 0 0 0 0 0 0 0 0 0 ...
## $ population : num 45195777 45195777 45195777 45195777 45195777 ...
## $ population_density : num 16.2 16.2 16.2 16.2 16.2 ...
## $ median_age : num 31.9 31.9 31.9 31.9 31.9 31.9 31.9 31.9 31.9 31.9 ...
## $ aged_65_older : num 11.2 11.2 11.2 11.2 11.2 ...
## $ aged_70_older : num 7.44 7.44 7.44 7.44 7.44 ...
## $ gdp_per_capita : num 18934 18934 18934 18934 18934 ...
## $ extreme_poverty : num 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 ...
## $ cardiovasc_death_rate : num 191 191 191 191 191 ...
## $ diabetes_prevalence : num 5.5 5.5 5.5 5.5 5.5 5.5 5.5 5.5 5.5 5.5 ...
## $ female_smokers : num 16.2 16.2 16.2 16.2 16.2 16.2 16.2 16.2 16.2 16.2 ...
## $ male_smokers : num 27.7 27.7 27.7 27.7 27.7 27.7 27.7 27.7 27.7 27.7 ...
## $ handwashing_facilities : num NA NA NA NA NA NA NA NA NA NA ...
## $ hospital_beds_per_thousand : num 5 5 5 5 5 5 5 5 5 5 ...
## $ life_expectancy : num 76.7 76.7 76.7 76.7 76.7 ...
## $ human_development_index : num 0.845 0.845 0.845 0.845 0.845 0.845 0.845 0.845 0.845 0.845 ...
Considerando el nivel de medición de las variables de nuestra base de datos, es necesario calcular los estadísticos de tendencia central que usted ya conoce, para todas las variables. Idealmente, se sugiere repetir el ejercicio para todas las variables.
# a) Media
mean(sud.covid$new_cases, na.rm = TRUE)
## [1] 5348.589
mean(sud.covid$new_cases_per_million, na.rm = TRUE)
## [1] 92.28003
#b) Mediana
median(sud.covid$new_cases, na.rm = TRUE)
## [1] 973
median(sud.covid$new_cases_per_million, na.rm = TRUE)
## [1] 67.036
# c) Moda: Como R no contiene una sintaxis específica para obtener la moda, es posible crear una función (objeto: función) para calcularla.
library(stats)
mode <- function(x) {
ux <- na.omit(unique(x) )
tab <- tabulate(match(x, ux)); ux[tab == max(tab) ]
}
# Luego, calcula la moda usando la función creada:
mode(sud.covid$new_cases_per_million)
## [1] 0
#Rango
min(sud.covid$new_cases_per_million, na.rm = TRUE) #na.rm permite remover o eliminar los casos perdidos (NA)
## [1] -450.772
max(sud.covid$new_cases_per_million, na.rm = TRUE)
## [1] 844.914
range(sud.covid$new_cases_per_million, na.rm = TRUE)
## [1] -450.772 844.914
#Varianza y Desviacion Estandar
var(sud.covid$new_cases_per_million, na.rm = TRUE)
## [1] 9931.873
sd(sud.covid$new_cases_per_million, na.rm = TRUE)
## [1] 99.65878
summary(sud.covid$new_cases_per_million) #Resumen estadístico
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -450.77 10.82 67.04 92.28 144.08 844.91 116
library(mosaic)
favstats(sud.covid$new_cases_per_million) #favstats entrega principales estadisticos para vectores.
## min Q1 median Q3 max mean sd n missing
## -450.772 10.82 67.036 144.0755 844.914 92.28003 99.65878 3887 116
En este caso, podemos hacer diversos cruces de variables, dado que disponemos de diversos vectores (más otros que se podrían crear usando estos mismos datos), lo que nos permite observar las relaciones que pudieran existir entre estas variables a través del coeficiente de correlación de Pearson. Para ello, primero utilizaremos la función plot() para observar la relación entre los distintos pares de variables. Luego podemos utilizar la función cor(), que nos permite calcular los coeficientes de Pearson para todos los posible cruces entre dos variables. Más información sobre estas dos funciones la pueden encontrar en (inglés) https://www.statmethods.net/stats/correlations.html y https://www.statmethods.net/graphs/scatterplot.html
¿Es posible observar una tendencia a partir de los Gráficos de Dispersión? (Completar para todos los pares de variables)
#Crea gráficos de dispersión entre dos vectores.
plot(sud.covid$new_cases_per_million, sud.covid$new_deaths_per_million, main="Gráfico de dispersión",
xlab="Casos diarios por millón hab.", ylab="Fallecimientos diarios por millón hab.", pch=19)
Considerando la cantidad de casos y algunos problemas de nuestra base de datos (p. ej.: días con casos negativos), necesitamos filtrar los outliers para poder observar de mejor forma las relaciones. Este procedimiento permite visualizar de mejor forma la relación entre las variables, pero debe utilizarse con mucha precaución, pues básicamente estaremos borrando aquellos reportes que incorporaban muchos nuevos casos y/o fallecimientos y que, por tanto, puede resultar ser información de gran valor. En este caso vamos a quedarnos con aquellos reportes diarios por país que tengan valores positivos en la cantidad de nuevos casos por millón, así como aquellos que tengan valores inferiores 600 en esta misma variable y menores a 30 en la cantidad diaria de fallecimientos.
sud.covid <- filter(sud.covid, new_cases_per_million >= 0 & new_cases_per_million <=600 & new_deaths_per_million <= 30)
plot(sud.covid$new_cases_per_million, sud.covid$new_deaths_per_million, main="Gráfico de dispersión",
xlab="Casos diarios por millón hab.", ylab="Fallecimientos diarios por millón hab.", pch=19)
¿Puedes observar alguna relación entre las variables? ¿Los puntos en el gráfico se encuentran agrupados o dispersos?
Como es probable que a veces los gráficos de dispersión no muestren con claridad una relación, el coeficiente de correlación resulta de gran ayuda, dado que permite cuantificar:
Para obtener el coeficiente de correlación para cada par de variables, se puede utilizar la función “cor()”, pero solo con vectores. Para ello vamos a analizar la relación entre 5 factores:
covid.cor <- subset(sud.covid, select = c(new_cases_per_million, new_deaths_per_million, human_development_index, gdp_per_capita, people_vaccinated_per_hundred))
cor(covid.cor, use="pairwise.complete.obs", method="pearson") #Pairwise nos va a permitir usar toda la información disponible para cada par de variables.
## new_cases_per_million new_deaths_per_million
## new_cases_per_million 1.0000000 0.76919469
## new_deaths_per_million 0.7691947 1.00000000
## human_development_index 0.2693649 0.15309612
## gdp_per_capita 0.1179776 -0.04684571
## people_vaccinated_per_hundred 0.3111324 0.06989306
## human_development_index gdp_per_capita
## new_cases_per_million 0.2693649 0.11797756
## new_deaths_per_million 0.1530961 -0.04684571
## human_development_index 1.0000000 0.78286701
## gdp_per_capita 0.7828670 1.00000000
## people_vaccinated_per_hundred 0.4089363 0.46973828
## people_vaccinated_per_hundred
## new_cases_per_million 0.31113238
## new_deaths_per_million 0.06989306
## human_development_index 0.40893634
## gdp_per_capita 0.46973828
## people_vaccinated_per_hundred 1.00000000
#Esta función la podemos correr automáticamente o guardarla dentro de un objeto como se describe a continuación.
correlacion <- cor(covid.cor, use="pairwise.complete.obs", method="pearson")
print(correlacion)
## new_cases_per_million new_deaths_per_million
## new_cases_per_million 1.0000000 0.76919469
## new_deaths_per_million 0.7691947 1.00000000
## human_development_index 0.2693649 0.15309612
## gdp_per_capita 0.1179776 -0.04684571
## people_vaccinated_per_hundred 0.3111324 0.06989306
## human_development_index gdp_per_capita
## new_cases_per_million 0.2693649 0.11797756
## new_deaths_per_million 0.1530961 -0.04684571
## human_development_index 1.0000000 0.78286701
## gdp_per_capita 0.7828670 1.00000000
## people_vaccinated_per_hundred 0.4089363 0.46973828
## people_vaccinated_per_hundred
## new_cases_per_million 0.31113238
## new_deaths_per_million 0.06989306
## human_development_index 0.40893634
## gdp_per_capita 0.46973828
## people_vaccinated_per_hundred 1.00000000
Finalmente, para aprovechar todas las potencialidades gráficas que ofrece R Studio, recomiendo utilizar una función transferida por el Prof. Juanjo Medina (U. Manchester) sobre el paquete pairs(), quien tiene su curso completo escrito en R Markdown y GitHub, abierto y disponible a todo público (pero en inglés): https://jjmedinaariza.github.io/R-for-Criminologists/
Pueden revisar la sintaxis de la función, para comprenderla y buscar más información en la web (hay mucha!). Es importante que primero corran toda la función antes de aplicarla. Básicamente crea las capas con las cuales reporta el coeficiente de Pearson (panel.cor) y los histogramas (panel.hist) para cada variable.
Una vez creada la función, pueden utilizarla para cualquier cruce entre vectores. Eventualmente podría utilizarse con variables ordinales tipo escala, pero en realidad eso no sería correcto, pues para ello se recomienda utilizar la correlación de Spearman o correlaciones policóricas.
#Funciones panel.cor y panel.hist
panel.cor <- function(x, y, digits=2, prefix="", cex.cor, ...){
usr <- par("usr")
on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y, use="pairwise.complete.obs", method="pearson")) #utilizaremos coef pearson
txt <- format(c(r, 0.123456789), digits=digits)[1]
txt <- paste0(prefix, txt, sep="")
if(missing(cex.cor)) cex.cor <- 0.9/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * (1+r)/2)
}
panel.hist <-function(x, ...){
usr <- par("usr")
on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5))
h <- hist(x, plot = FALSE)
breaks <- h$breaks
nB <- length(breaks)
y <- h$counts
y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y, col="lightskyblue1", ...)
}
Ahora vamos a utilizar esas funciones en el gráfico que nos ofrece pairs():
#Aplicación
pairs(covid.cor, upper.panel=panel.smooth, lower.panel=panel.cor, diag.panel=panel.hist)
A partir de lo anterior, es posible concluir respecto a la importancia de reducir la cantidad de nuevos contagios para un efectivo control de la pandemia y una disminución de la tasa de letalidad. A la vez, resulta necesario considerar una estrategia regional de inmunización que no condicione el acceso al proceso de vacunación al nivel de riqueza o desarrollo de cada país.
Utilizando la misma base de datos, para esta semana deberá revisar la correlación que existe entre la cantidad de personas vacunadas (por cien) y otras 4 variables que no hayamos utilizado previamente (a elección). Para ello deberá realizar el análisis correlacional siguiendo los pasos descritos anteriormente y responder las consultas que se plantean en la siguiente sección.
La UNAP cuenta con algunas estaciones de monitoreo de las condiciones metereológicas en diversos puntos de la región. En el caso de Iquique y Alto Hospicio, la estación se encuentra en el Campus Huayquique de la UNAP, tal como se muestra en la siguiente ficha informativa de la estación.
Como se puede observar a continuación, la base de datos original contiene información a nivel de horas de cada día. Es decir, cada fila está representando a una hora de cada uno de los días entre el 01/ENE/2020 y el 15/JUN/2020. Por lo mismo vamos a manipular levemente los datos, para que puedanla interpretación sea más simple, a nivel de día.
#Importar los datos desde Excel
library(readxl)
Huayquique2020 <- read_excel("~/Google Drive/UNAP/Cursos/Estadistica/Meteorología UNAP/Huayquique2020.xlsx",
col_types = c("date", "date", "numeric",
"numeric", "numeric", "numeric"))
#Observemos los primeros casos, para hacernos una idea de los datos.
head(Huayquique2020)
## # A tibble: 6 x 6
## Date Time Temp Hum Speed Bar
## <dttm> <dttm> <dbl> <dbl> <dbl> <dbl>
## 1 2020-01-01 00:00:00 1899-12-31 00:00:00 21.2 75 0 1010.
## 2 2020-01-01 00:00:00 1899-12-31 01:00:00 20.8 76 0 1011.
## 3 2020-01-01 00:00:00 1899-12-31 02:00:00 20.4 78 0 1010.
## 4 2020-01-01 00:00:00 1899-12-31 03:00:00 20.3 78 0 1010.
## 5 2020-01-01 00:00:00 1899-12-31 04:00:00 20 79 0 1009.
## 6 2020-01-01 00:00:00 1899-12-31 05:00:00 20.1 79 0 1008.
#Podemos extraer la hora de una variable de fecha, con el paquete 'lubridate'
library(lubridate)
Huayquique2020$hora <- hour(Huayquique2020$Time)
#Creamos un data.frame que contiene la media diaria para cada variable con 'dplyr'
library(dplyr)
hqq <- Huayquique2020 %>%
group_by(Date) %>%
summarise(temp=mean(Temp),hum=mean(Hum),viento=mean(Speed),pres=mean(Bar))
#Reducimos los decimales para reducir la complejidad (no siempre es aconsejable hacerlo, depende del rango de la variable)
hqq$temp <- round(hqq$temp, digits = 1)
hqq$hum <- round(hqq$hum, digits = 1)
hqq$viento <- round(hqq$viento, digits = 2)
hqq$pres <- round(hqq$pres, digits = 1)
#Exportamos a CSV, para subirlo a GitHub
write.csv2(hqq,"/Users/danielquinterosr/Google Drive/UNAP/Cursos/Estadistica/Clases R/huayquique.csv", row.names = FALSE, fileEncoding = "UTF-8") #UTF-8 permite reconocer acentos y caracteres específicos del castellano, como la Ñ
Ahora ya contamos con los datos necesarios y en el formato adecuado para poder trabajar. Esta preparación del conjunto de datos que se va a utilizar posteriormente, es un paso muy importante en el proceso de análisis estadístico. Importaremos los datos desde GitHub:
hqq <- read.csv2(url("https://raw.githubusercontent.com/danielquinterosr/IntroR/master/huayquique.csv"), header=TRUE, sep=";", fileEncoding = "UTF-8")
str(hqq)
## 'data.frame': 167 obs. of 5 variables:
## $ Date : chr "2020-01-01" "2020-01-02" "2020-01-03" "2020-01-04" ...
## $ temp : num 21.8 21.9 21.3 21.2 22.2 22.7 22.3 21.6 20.9 20.5 ...
## $ hum : num 69.7 69.2 74.2 76.6 72.4 68.3 74 72 72.9 75.2 ...
## $ viento: num 0.29 0.63 0.71 0.47 0.32 0.31 0.44 0.81 0.96 0.5 ...
## $ pres : num 1009 1009 1011 1011 1011 ...
Ahora que ya tenemos cargados los datos en nuestro entorno de R, deberás aplicar los estadísticos que hemos revisado en esta unidad a cada una de las variables contenidas en el data.frame hqq (temperatura, humedad, viento y presión atmosférica -bars-). Sugiero siempre comenzar con una mirada general al conjunto de datos, utilizando la función str().
Recuerda que el primer paso siempre consistirá en realizar una exploración descriptiva de las variables que vayamos a utilizar en nuestro análisis.
library(mosaic)
favstats(hqq$temp)
## min Q1 median Q3 max mean sd n missing
## 15.6 18.8 20.9 21.9 24.1 20.39042 1.923139 167 0
favstats(hqq$hum)
## min Q1 median Q3 max mean sd n missing
## 67.7 72.25 73.8 76 84 74.17904 3.120566 167 0
favstats(hqq$viento)
## min Q1 median Q3 max mean sd n missing
## 0 0.23 0.41 0.545 1.01 0.3986228 0.2198053 167 0
favstats(hqq$pres)
## min Q1 median Q3 max mean sd n missing
## 1004.6 1008.4 1009.2 1010.95 1014.3 1009.425 2.059293 167 0
Luego podemos comenzar a explorar las relaciones entre dos variables, para lo cual tomaremos temperatura y viento ¿Puedes observar alguna relación?
plot(hqq$temp, hqq$viento, main="Gráfico de dispersión",
xlab="temperatura", ylab="viento", pch=19)
El ojo humano tiene una serie de limitaciones para poder observar gráficamente las relaciones. Por fortuna, tanto el entrenamiento en análisis de datos como las herramientas de apoyo nos pueden ser de gran utilidad. En este caso, podemos agregar al gráfico de dispersión el trazado de la curva de regresión para observar de mejor forma la posible relación. Esta curva nos será muy valiosa más adelante cuando veamos regresiones lineales, pero por mientras puedes ver el texto de Ferris Ricthey (Cap.14) y aprovechar la información que entrega:
plot(hqq$temp, hqq$viento, main="Gráfico de dispersión",
xlab="temperatura", ylab="viento", pch=19)
abline(lm(hqq$viento~hqq$temp), col="red") # regression line (y~x)
Veamos ahora la relación entre temperatura y humedad:
plot(hqq$temp, hqq$hum, main="Gráfico de dispersión",
xlab="temperatura", ylab="humedad", pch=19)
abline(lm(hqq$hum~hqq$temp), col="red") # regression line (y~x)
Y la relación entre temperatura y presión atmosférica:
plot(hqq$temp, hqq$pres, main="Gráfico de dispersión",
xlab="temperatura", ylab="presión", pch=19)
abline(lm(hqq$pres~hqq$temp), col="red") # regression line (y~x)
hqq.cor <- subset(hqq, select = c(temp,hum,viento,pres))
cor(hqq.cor, use="pairwise.complete.obs", method="pearson") #Pairwise nos va a permitir usar toda la información disponible para cada par de variables.
## temp hum viento pres
## temp 1.0000000 -0.25101423 0.29913839 -0.6057396
## hum -0.2510142 1.00000000 -0.04097778 0.1778667
## viento 0.2991384 -0.04097778 1.00000000 -0.1996467
## pres -0.6057396 0.17786673 -0.19964666 1.0000000
pairs(hqq.cor, upper.panel=panel.smooth, lower.panel=panel.cor, diag.panel=panel.hist)
Hasta acá hemos trabajado sólo con datos poblacionales. En el caso de los datos de la pandemia, cada fila representaba un reporte diario en una locación determinada. Por su parte, los datos meteorológicos fueron resumidos a partir de sus valores diarios, con lo cual usando la función cor() es suficiente para determinar la existencia, dirección y fuerza de la relación entre dos variables.
Sin embargo, como recordarás del curso sobre estadística inferencial, la gran diferencia entre trabajar con datos poblacionales y datos muestrales, es que en este último caso, tenemos la posibilidad cierta de cometer un error de estimación. ¿Cómo podemos asegurarnos, entonces, de que las relaciones que observemos en nuestras muestras sean efectivamente los datos probables de la población?
J. Bruna & F. Reyes (2020) https://danielquinterosr.github.io/EST2020.html
Según la Organización Mundial de la Salud, la depresión es comprendida como un trastorno mental, que se caracteriza por la presencia de tristeza, pérdida del interés o placer, sentimientos de culpa o falta de autoestima, sensación de cansancio y falta de concentración. En su forma más grave puede conducir al suicidio. En esta línea, la Salud Mental es comprendida como un estado de completo bienestar físico, mental y social, no solamente la ausencia de afecciones o enfermedades, sino que también el detectarlas, tratarlas y finalmente rehabilitarse de aquella enfermedad para el conseguimiento de un desarrollo pleno.
En este sentido, el presente trabajo buscó conocer, identificar y analizar la intensidad e impacto de los problemas psicológicos que se presentan en la población de la Región de Tarapacá en tiempos de pandemia por COVID-19.
A partir de las diferencias que muestra el gráfico, podemos observar que las mujeres son las que experimentan con mayor frecuencia estas emociones negativas, y estadísticamente son las mujeres las que tienen mayor carga laboral, domestica, familiar, entre otras. Es por esto que creemos que son las mujeres las más afectadas en salud mental debido a la pandemia COVID-19.
est2020 <- read.csv(file = ("https://raw.githubusercontent.com/danielquinterosr/EncuestaSocialTarapaca2020/master/est2020.csv"), header=TRUE, sep=";", fileEncoding = "UTF-8")
library(dplyr)
est.s <- est2020 %>%
filter(sexo!="Otro")
est.s$sexo <- factor(est.s$sexo,levels = c("Mujer","Hombre"))
est.s <-subset(est2020, select = c (d1.angustia, d1.ansiedad, d1.estres, d1.frustracion, d1.incertidumbre, d1.enojo, d1.miedo, d1.soledad, d1.tristeza, d2.atencion, d3.farmacos, d4.tratamiento))
est2020$d1.angustia <- factor(est2020$d1.angustia,
levels = c("Nunca",
"Casi nunca",
"A veces",
"Frecuentemente",
"Siempre"))
est2020$d1.ansiedad <- factor(est2020$d1.ansiedad,
levels = c("Nunca",
"Casi nunca",
"A veces",
"Frecuentemente",
"Siempre"))
est2020$d1.estres <- factor(est2020$d1.estres,
levels = c("Nunca",
"Casi nunca",
"A veces",
"Frecuentemente",
"Siempre"))
est2020$d1.frustracion <- factor(est2020$d1.frustracion,
levels = c("Nunca",
"Casi nunca",
"A veces",
"Frecuentemente",
"Siempre"))
est2020$d1.incertidumbre <- factor(est2020$d1.incertidumbre,
levels = c("Nunca",
"Casi nunca",
"A veces",
"Frecuentemente",
"Siempre"))
est2020$d1.enojo <- factor(est2020$d1.enojo,
levels = c("Nunca",
"Casi nunca",
"A veces",
"Frecuentemente",
"Siempre"))
est2020$d1.miedo <- factor(est2020$d1.miedo,
levels = c("Nunca",
"Casi nunca",
"A veces",
"Frecuentemente",
"Siempre"))
est2020$d1.soledad <- factor(est2020$d1.soledad,
levels = c("Nunca",
"Casi nunca",
"A veces",
"Frecuentemente",
"Siempre"))
est2020$d1.tristeza <- factor(est2020$d1.tristeza,
levels = c("Nunca",
"Casi nunca",
"A veces",
"Frecuentemente",
"Siempre"))
est2020$d2.atencion <- factor(est2020$d2.atencion,
levels = c("No",
"Tal vez",
"Si"))
est2020$d3.farmacos <- factor(est2020$d3.farmacos,
levels = c("No",
"Si"))
est2020$d4.tratamiento <- factor(est2020$d4.tratamiento,
levels = c("No estaba en ningún tratamiento antes de la pandemia",
"Si, con total normalidad",
"Si, pero con algunas dificultades",
"No he podido continuar el tratamiento"))
library(dplyr)
est.s <- est2020 %>%
filter(sexo!="Otro")
est.s$sexo <- factor(est.s$sexo,levels = c("Mujer","Hombre"))
est.s$sexo <- as.factor(est.s$sexo)
est.s1 <- subset(est.s, select = c (d1.angustia,
d1.ansiedad,
d1.estres,
d1.frustracion,
d1.incertidumbre,
d1.enojo,
d1.miedo,
d1.soledad,
d1.tristeza))
names(est.s1)= c('"Ha sentido angustia"',
'"Ha sentido ansiedad"',
'"Ha sentido estres"',
'"Ha sentido frustracion"',
'"Ha sentido incertidumbre"',
'"Ha sentido enojo"',
'"Ha sentido miedo"',
'"Ha sentido soledad"',
'"Ha sentido tristeza"')
library(likert)
## Loading required package: xtable
##
## Attaching package: 'likert'
## The following object is masked from 'package:dplyr':
##
## recode
likert.s1 <- likert(est.s1, grouping=factor(est.s$sexo,levels = c("Mujer","Hombre")))
likert.s1g <- likert(summary = likert.s1$results, grouping = likert.s1$results[,1])
likert.s1g <- likert(summary = likert.s1g$results[-c(2)])
likert.s1g <- likert(summary = subset(likert.s1g$results, select = -c(Group)),grouping = likert.s1$results[,1])
library(rcartocolor)
paletalikert2 = carto_pal(5, "TealGrn")
graf.smental <- plot(likert.s1g,center=3,plot.percent.neutral = TRUE,include.center = TRUE,plot.percents=TRUE,text.color = "black",plot.percent.low = TRUE, plot.percent.high = TRUE,text.size=5) +
scale_fill_manual(values=paletalikert2, breaks=likert.s1g$levels, labels=likert.s1g$levels, guide = guide_legend(reverse = TRUE)) +
theme(plot.title = element_text(size = 22, face = "bold", hjust = 0),
legend.title = element_blank(),
legend.key.size = unit(0.9,"line"),
legend.position="right",
legend.text = element_text(size = 12),
plot.caption = element_text(size=14),
plot.subtitle=element_text(size=14, face = "italic",color = "#6b6e6c"),
strip.text=element_text(size=13),
axis.text = element_text(size=12),
strip.background = element_rect(size=12),
rect = element_rect(fill = "transparent"),
panel.background = element_rect(colour = 'transparent'),
plot.background = element_rect(fill = "transparent", color = NA)) +
labs(title="¿Con qué frecuencia diría usted que ha sentido \n las siguientes emociones durante la actual pandemia del Covid-19?",y="%",caption="(N=244)",subtitle = '\nENCUESTA SOCIAL TARAPACÁ 2020 \n \n (Nunca) (A veces + Frecuentemente + Casi siempre + Siempre)')
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
#ggsave(graf.smental, filename = "likert-saludmental.png",width = 30, height = 40, units = "cm") guarda como png (recomendado)
#print(graf.smental) imprime directamente (no recomendado)
La matriz con los coeficientes de correlación y los gráficos de dispersión dieron cuenta de una fuerte correlación entre todas las variables. Según las pruebas realizadas usando el coeficiente de correlación de Spearman, las siguientes fueron estadísticamente significativas:
Miedo – Angustia (0.74): Esta correlación monótona va en aumento, es decir, cuando la variable “x” aumenta y la variable “y” nunca disminuye. Se observa una fuerte correlación entre ambos, a mayor miedo, mayor angustia, a mayor angustia, mayor miedo.
Angustia – Ansiedad (0.63): En esta correlación podemos ver que existe una fuerte asociación de rango, que da cuenta de una correlación entre ambas emociones.
Soledad – Tristeza (0.60): En esta correlación existe la misma dinámica de las anteriores, en donde se ve que la variable en soledad mientras va aumentando, la triste de igual manera se encuentra creciendo. Y viceversa.
Tristeza – Miedo (0.58): Estas dos emociones triste y miedo se correlación pero su correlación no es tan fuerte como la de las anteriores.
Ansiedad – Soledad (0.50): Esta correlación muestra la relación entre ansiedad y soledad *(puedo sentirme solo/a y sentir ansiedad o viceversa).
#Funciones panel.cor y panel.hist
panel.cor <- function(x, y, digits=2, prefix="", cex.cor, ...){
usr <- par("usr")
on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y, use="pairwise.complete.obs", method="spearman")) #utiliza coef rho spearman
txt <- format(c(r, 0.123456789), digits=digits)[1]
txt <- paste0(prefix, txt, sep="")
if(missing(cex.cor)) cex.cor <- 0.9/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * (1+r)/2)
}
panel.hist <-function(x, ...){
usr <- par("usr")
on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5))
h <- hist(x, plot = FALSE)
breaks <- h$breaks
nB <- length(breaks)
y <- h$counts
y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y, col="lightskyblue1", ...)
}
est.s <- est2020 %>%
filter(sexo!="Otro")
est.s$sexo <- factor(est.s$sexo,levels = c("Mujer","Hombre"))
est.s$sexo <-factor(est.s$sexo)
est.s5 <- subset(est.s, select = c (d1.angustia,
d1.ansiedad,
d1.estres,
d1.frustracion,
d1.incertidumbre,
d1.enojo,
d1.miedo,
d1.soledad,
d1.tristeza))
est.s5$d1.angustia <-as.numeric(est.s5$d1.angustia)
est.s5$d1.ansiedad <-as.numeric(est.s5$d1.ansiedad)
est.s5$d1.estres <-as.numeric(est.s5$d1.estres)
est.s5$d1.frustracion <-as.numeric(est.s5$d1.frustracion)
est.s5$d1.incertidumbre <-as.numeric(est.s5$d1.incertidumbre)
est.s5$d1.enojo <-as.numeric(est.s5$d1.enojo)
est.s5$d1.miedo <-as.numeric(est.s5$d1.miedo)
est.s5$d1.soledad <-as.numeric(est.s5$d1.soledad)
est.s5$d1.tristeza <-as.numeric(est.s5$d1.tristeza)
est.s5c <- subset(est.s5, select = c (d1.angustia,
d1.ansiedad,
d1.enojo,
d1.miedo,
d1.soledad,
d1.tristeza))
names(est.s5c)=c("Angustia",
"Ansiedad",
"Enojo",
"Miedo",
"Soledad",
"Tristeza")
cor(est.s5c, use="pairwise.complete.obs", method="spearman")
## Angustia Ansiedad Enojo Miedo Soledad Tristeza
## Angustia 1.0000000 0.6373401 0.5998351 0.7490110 0.5547634 0.6385509
## Ansiedad 0.6373401 1.0000000 0.5081071 0.5141677 0.5037798 0.5175637
## Enojo 0.5998351 0.5081071 1.0000000 0.4969305 0.5180695 0.5585789
## Miedo 0.7490110 0.5141677 0.4969305 1.0000000 0.5318824 0.5825093
## Soledad 0.5547634 0.5037798 0.5180695 0.5318824 1.0000000 0.6046055
## Tristeza 0.6385509 0.5175637 0.5585789 0.5825093 0.6046055 1.0000000
Como conclusión final creemos que se les debe de poner atención a la salud mental de las personas, que se generen programas y oportunidades para que las personas de la región, especialmente las mujeres, puedan mejorar su salud mental. Creemos también, por los resultados, que la población no reconoce que necesitan atención, pero que aun así han sentido y pasado por emociones negativas. Queremos que generen espacios y ambientes agradables en su hogar, que puedan organizarse y planificarse, que tengan tiempo libre y realicen algo que les guste. Queremos que la salud mental de los tarapaqueños mejore y se tome en cuenta.
pairs(est.s5c, pch=".", upper.panel=panel.smooth, lower.panel=panel.cor, diag.panel=panel.hist, main = "Salud Mental - Coeficientes de correlación (Rho)")
Cuando queremos realizar un proceso de inferencia estadística, por ejemplo, para observar si la relación observada en la muestra entre miedo y angustia puede ser replicada a nivel poblacional, siempre deberemos realizar una prueba de hipótesis que, en este caso, adopta la siguiente forma:
cor.test(est.s5c$Angustia, est.s5c$Miedo, use="pairwise.complete.obs", method="spearman")
## Warning in cor.test.default(x, y, ...): Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: x and y
## S = 607667, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.749011
#para mirar otras posibilidades de la sintaxis, puedes utilizar el siguiente comando para abrir el archivo de ayuda:
#?cor.test
Como puedes observar, el valor de significancia (p-value) es menor al puntaje crítico de 0,05 que se corresponde con un nivel de confianza de 95% (podrías utilizar otro puntaje), con lo cual existe evidencia estadísticamente significativa como para rechazar la hipótesis nula y establecer que las variables miedo y angustia presentan una fuerte asociación en la población o, en palabras de Hernández Sampieri et al. (2010), una correlación positiva considerable. Pero recuerda: correlación no es causalidad.
© danielquinterosr